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Abstract 

The fractal structure of high-temperature graphs of the three-dimensional Ising and XY models 
is investigated by simulating these graphs directly on a cubic lattice and analyzing them with the 
help of percolation observables. The Ising graphs are shown to percolate right at the Curie critical 
point. The diverging length scale relevant to the graphs in the vicinity of the percolation threshold is 
shown to be provided by the spin correlation length. The fractal dimension of the high-temperature 

£ 

graphs at criticality is estimated to be D = 1.7349(65) for the Ising and D = 1.7626(66) for the 
XY model. 
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I. INTRODUCTION 



The high-temperature (HT) expansion is a powerful tool to study the critical properties 
of lattice spin models pp. In this approach, the partition function and correlation functions 
are calculated by counting graphs on the lattice with each graph representing a certain 
contribution. Traditionally, such an expansion is carried out exactly to a given order by 
enumerating all possible ways a graph of given size and topology can be drawn on the 
lattice. This exact approach, involving combinatorial and graph-theoretical algorithms, is 
notoriously challenging and laborious, with each additional order requiring typically about 
the same amount of effort needed for all previous orders combined. 

We have developed a different approach |2J in which the HT representation of lattice spin 
models is studied by means of Monte Carlo simulations. The HT graphs along the links of 
the underlying lattice are generated through a Metropolis plaquette update that proposes 
a local change in the existing graph configuration. At high temperatures, only a few small 
graphs generated this way can be found scattered throughout the lattice. As the temperature 
is lowered, graphs start to fill the lattice by growing larger and becoming more abundant. At 
temperatures below the critical temperature, the lattice becomes filled with graphs. A typical 
graph configuration now consists of one big graph spanning the entire lattice and a collection 
of much smaller graphs (see Fig. [Tj). The steady increase in the number of occupied links and 
the appearance of graphs spanning the lattice as the temperature is lowered are reminiscent of 
a percolation process. The use of percolation observables therefore suggests itself to analyze 
the graph configurations. For these observables to have bearing on the critical properties 
of the model under investigation, it is necessary that the HT graphs percolate right at the 
thermal critical point. For the Ising model on a square lattice we numerically showed that the 
percolation threshold indeed coincides with the (exactly known) Curie temperature [2]. In 
other words, the phase transition in this lattice spin model to the ordered, low-temperature 
state manifests itself through a proliferation of HT graphs. Moreover, the fractal structure 
of closed and open graphs was shown to encode the standard critical exponents [2j [3j HI [5] . 

The purpose of this paper is to extend this geometric Monte Carlo approach to three 
dimensions. Two-dimensional (2D) spin models arguably form a very special class of models, 
in particular the Ising model as it is self-dual. It therefore is a priori not obvious that this 
geometric Monte Carlo approach together with the use of percolation observables is viable 
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FIG. 1: Distribution of Ising HT graphs on a cubic lattice of linear size L = 24 at the percolation 
threshold. Note the presence of a single big graph and many much smaller graphs. 

in 3D. 

Another Monte Carlo algorithm for studying HT representations of classical statistical 
models has been put forward by Prokof 'ev and Svistunov [6 J . That so-called worm algorithm 
is much more efficient than the conventional local update we use. The dynamic exponent 
z characterizing the divergence of the autocorrelation time r when the critical point is 
approached, r ~ £ 2 , with £ the correlation length, is close to zero for the worm algorithm 
while it is larger than two for the plaquette update we use. The plaquette update, on the 
other hand, has the virtue that it provides a direct and clean implementation of the HT 
representations of the spin models we consider. Prokof 'ev and Svistunov [I] recently applied 
their algorithm to the 3D complex |0| 4 theory to determine the fractal dimension of the 
HT graphs in that theory. These graphs are unlike those in the XY model as links and 
vertices carry different weights in the two models. Moreover, because the update algorithms 
differ, the graphs in the two models are simulated in completely different ways. Nevertheless, 
since the |0| 4 theory is in the same universality class as the XY model, both types of HT 
graphs should yield the same fractal dimension. We set out below to investigate whether 
universality holds for the fractal structure of HT graphs. 

The paper is organized as follows. The next section introduces the Metropolis plaquette 
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update algorithm used in this Monte Carlo study together with the percolation observables 
applied to analyze the HT graphs. The subsequent two sections present our results for the 
Ising (Sec. Ill) and XY (Sec. IV) models, and the paper ends in Sec. [Vjwith a summary and 
conclusions. 



II. SIMULATION AND DATA ANALYSIS TECHNIQUES 

To be specific, we consider the HT representation of 0(N) lattice spin models described 
by the Hamiltonian 

# = -J^S x -S x ,, (1) 

<x,x'> 

with the interaction, characterized by the parameter J, restricted to spins on nearest- 
neighbor sites, so that the sum in Eq. Q runs only over nearest-neighbor pairs. The spin 
variable S x = (S 1 *, ■ ■ ■ , S$?) located at each site x of the cubic lattice has a fixed length, 
S x = 1. Simulations are carried out for the 3D Ising (N = 1) and XY (N = 2) models. 



A. Ising Model 

The HT representation of the 3D Ising model on a cubic lattice with periodic boundary 
conditions consisting of N sites and 3 N links [7] , 

Z = (cosh/3) 3N 2 N Rb > ( 2 ) 

closed 
graphs 

provides an alternative, but completely equivalent description of the spin model. In this 
representation, which is purely geometric in nature, spin degrees of freedom are swapped for 
link variables. The representation ^ of the partition function can be visualized as a sum 
over all possible closed graphs that can be drawn on the lattice. Each occupied link carries 
a factor K = tanh/3, with (3 the inverse temperature, where for convenience the coupling 
constant J is set to unity. In the entire temperature range < (3 < oo, < K < 1. The 
minimum number of occupied links h needed to form a closed graph is four on a cubic lattice. 
The internal energy 

is determined by the average number (b) of occupied links. 



The central idea of our geometric Monte Carlo approach |2J is to directly simulate the 
graphs contributing to the partition function. The HT representation (p| suggests the fol- 
lowing local Metropolis update algorithm [8]. 

The probability distribution P(G) for a given graph configuration G reads in equilibrium 

P{G) = i(cosh/5) 37V 2 7V K b . (4) 

Zj 

Such a configuration can be reached from the configuration present after, say, t iterations in 
the following way 

P t+1 (G) = P t (G) + [Pt{G')W{G' - G) - P t (G)W(G - G')} , (5) 

G> 

where W(G — > G') is the probability for the system to move from the graph configuration G 
with b occupied links to the graph configuration G' with b' occupied links. In equilibrium, 
P t+ \{G) = Pt{G) = P(G), and the system satisfies detailed balance 

P(G')W(G' -> G) = P{G)W{G -> G'), (6) 

or 

W(G -> G") _ P(G') _ 

W{G' ^G)~ P{G) ~ ~K*' ^' 
As is custom with Metropolis algorithms, the acceptance rate pht of a proposed update 
is maximized by giving the largest of the two transition probabilities W(G — > G 1 ) and 
W(G' — > G) appearing in the ratio the largest possible value, which is one. That is, if the 
number of links b' in the proposed configuration is larger than the number of links b in the 
existing configuration, so that K b '/K b < 1, W(G' -> G) = 1 and W(G -> G') = p HT . If, on 
the other hand, b' < b, the proposed configuration carries a larger weight than the existing 
one and will always be accepted. 

The HT graphs are generated by taking the smallest possible closed graphs on the lattice, 
i.e., plaquettes, as building blocks. During a sweep of the lattice, all plaquettes are visited 
in a regular, typewriter fashion. For the Ising model, reflecting the underlying Z 2 spin 
symmetry, a link can either be empty or occupied. The links of a plaquette considered for 
update are changed from empty to occupied and vice versa (see Fig. [2] for an illustration). 
This is easily implemented by means of the binary rules + 1 = 1, 1 + 1 = 0, respectively. 
The acceptance rate of a proposed update reads in formula [8] 
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FIG. 2: An existing Ising HT graph on a cubic lattice (left panel) is updated with the help of a 
chosen plaquette (middle panel) into a new graph (right panel) . 

f K b '- b if V > b 
Put = < (8) 
I 1 else. 

The resulting and existing number of occupied links, b' and b, respectively, are related via 

b' = b + 4 - 2b n , (9) 

with bu denoting the number of links on the plaquette already occupied. By taking plaquettes 
as building blocks, the resulting graphs are automatically closed. 

Table [I] gives a summary of the number Nmc of Monte Carlo sweeps of the lattice of 
size L used for data collection in the temperature interval [/5i,/5 2 ], with N c &\ sweeps used 
for equilibration. The temperature intervals are sampled with i equidistant points. The 
largest lattice taken for temperature-dependent runs is L = 48, while runs at the percolation 
threshold are carried out on lattices up to L = 64. Since the plaquette update is a local 
update, autocorrelation times grow very large on large lattices. Analyzing the time series 
of the percolation strength, we estimate the autocorrelation time to vary from r »s 25 for 
L = 10 to as long as r ~ 2500 for L = 64. To somewhat reduce correlations between 
successive data points, we take measurements every fifth sweep of the lattice. Statistical 
errors are estimated by means of jackknife binning. Fits are carried out by using the nonlinear 
Marquardt-Levenberg algorithm for minimization of error weighted least squares. 
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TABLE I: Overview of parameters used in the simulations of the 3D Ising model. The temperature 
intervals \fli, P2} are sampled with i equidistant points. For each sampling point, A^mc Monte Carlo 
sweeps of the lattice of size L are used for data collection, with iV ca i sweeps used for equilibration. 

01 02 i L N MC iV ca i 

0.218000 0.233000 70 12 460000 60000 

0.170000 0.218000 190 12 120000 20000 

0.233000 0.270000 120 12 120000 20000 

0.220000 0.230000 60 16 460000 60000 

0.170000 0.220000 190 16 120000 20000 

0.230000 0.270000 140 16 120000 20000 

0.220500 0.229000 50 20 460000 60000 

0.170000 0.220500 190 20 120000 20000 

0.229000 0.270000 150 20 120000 20000 

0.220500 0.227100 40 24 460000 60000 

0.170000 0.220500 190 24 120000 20000 

0.227100 0.270000 160 24 120000 20000 

0.220500 0.226000 35 28 460000 80000 

0.170000 0.220500 180 28 120000 30000 

0.226000 0.270000 165 28 120000 30000 

0.220500 0.224900 30 32 460000 80000 

0.170000 0.220500 190 32 120000 30000 

0.224900 0.270000 170 32 120000 30000 

0.220500 0.223800 25 40 500000 120000 

0.170000 0.220500 180 40 130000 40000 

0.223800 0.270000 175 40 130000 40000 

0.220500 0.223500 25 48 560000 240000 

0.170000 0.220500 150 48 140000 50000 

0.223500 0.270000 160 48 140000 50000 
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B. XY Model 



For the XY model (N = 2), where the spins take values along a circle, the general spin 
Hamiltonian (JIJ with J = 1 reduces to 



H 



£ cos ( 

(x,x'> 



(10) 



with 6 X being the planar angle of the spin at site x relative to a fixed, but arbitrary axis. 
The HT expansion is facilitated by the use of the Fourier expansion 



3 /3cos(0) 



£ 



(ii) 



l=— oo 



where Ii(x) is the modified Bessel function of the first kind. The partition function, 



n 



2tT 



n 



-/3cos(0 x -6> x ,) 



(12) 



J (x,x'> 



then takes the well-known form [9] 



" d^x 



II E J ^C9)^ ( 



(13) 



involving the integers l XjX > defined on the links connecting the nearest neighbor sites x and 
x'. The spin degrees of freedom are now easily integrated out with the result 



(14) 



Here, the prime on the sums is to indicate that only configurations satisfying the zero 
divergence condition, J^ x , Z XjX / = at each site x contribute, where the sum runs over 
all nearest neighbors of x. It is generally accepted that the link variables i xx / at the right 



hand of Eq. (14) can be restricted to the values ±1 without changing the universality class. 



The partition function of the XY model can then be cast in a form analogous to the HT 
representation ^ of the Ising model, 



Z = I 3 N (P) £ K b 



(15) 



closed 
oriented graphs 



where K = Ix((3) / Io((3) , with < K < 1 for all /3, and use is made of the property that 
= I\{x). Since the link variable of the truncated model can take two nontrivial 
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values ±1, the graphs now have, in contrast to the Ising model, an orientation. A plaquette 
considered for update must therefore also be given a (randomly chosen) orientation. As 
binary rules we now have —1 + 1 = 0, + 1 = 1, — 1 = —1, 1 — 1 = in addition to 
the restrictions —1 — 1= 0, 1 + 1 = of singly occupancy which, as already mentioned, we 
expect not to change the universality class of the model. Apart from these modifications, 
the HT representation of the XY model can be handled in the same way as that of the Ising 
model. 

C. Observables 

The HT graphs are analyzed with the help of standard percolation observables [10J. An 
important characteristic is whether a graph spans the lattice or not. We say a graph does 
so already if it spans the lattice in just one of the three possible directions. By recording 
this each time the graphs are analyzed, one obtains the percolation probability Ps, which 
tends to zero in the limit (3 — > and to unity in the opposite limit (3 — > oo. Another 
important observable is the graph distribution n^, which gives the average number of graphs 
of b occupied links normalized by the volume. Close to the percolation threshold it assumes 
the form 

n b ~b- T e~ e \ 6 oc \K/K PCI - 1| 1/<T . (16) 

The exponents r and a are related to the fractal dimension D of the HT graphs via r = 
d/D + 1 with d = 3 the dimension of the space box, and [3 J 

a = l/uD, (17) 

where v is the exponent characterizing the divergence of the correlation length £ as the 
percolation threshold is approached, £ oc \K/K per — l\~ v '. An additional observable we 
measure is the percolation strength P^, which is defined as the size of the largest graph 
normalized by the volume. Finally, we also record the average graph size xg- 

The percolation threshold and the fractal dimension of the HT graphs follow from ap- 
plying finite-size scaling to these observables. According to scaling theory, the percolation 
probability and strength for different values of the tuning parameter K and for different lat- 
tice sizes L do not depend on these variables separately, but depend on them in a convoluted 
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way [TT] 

P S (*Q = Ps [L^(K/K pei - 1)] , P^K) = L~^' v Poo [L 1 ^(K/K PCI - 1)] . (18) 

Here, (5q determines the scaling dimension of the percolation strength, which is related to 
the fractal dimension D of the HT graphs through [10] 



D = d - p G /v. (19) 
The scaling dimension of the percolation probability is zero. The first scaling relation in 



Eq. (18) implies that the curves Ps(K) measured on lattices of different size all cross at the 
same point. This point, being volume independent, marks the percolation threshold K peT on 
the infinite lattice. That scaling relation in addition implies that the curves collapse onto 
a universal curve when replotted as a function of l}l v {Kj K pcr — 1). Similarly, the second 



scaling relation in Eq. ( 18 ) implies that if in addition to the horizontal axis also the vertical 
axis is properly rescaled, with the correct value for the ratio (3q/v, also the P^K) data fall 
onto a universal curve. That scaling relation in addition implies that measurements at the 
percolation threshold scale as P^K^) oc L _/3g//!/ , providing a good means of determining 



III. ISING MODEL RESULTS 

We start by simulating the HT graphs of the Ising model on a cubic lattice. Figure [3] 
shows the internal energy as obtained through Eq. (|3]) by measuring the density of occupied 
links. Apart from a small interval around the critical temperature, E is seen to be almost 
independent of volume, implying that the correlation length is much smaller here than the 
linear size of the smallest lattice considered (L = 12). When entering the critical region, 
the correlation length becomes larger and eventually exceeds the size of the largest lattice 
considered [L = 48). 

We first determine the location of the percolation threshold. To this end, we measure 
the percolation probability on lattices of different linear size L as a function of the inverse 
temperature (see Fig. [4]). The figure shows that the curves connecting the data points 
all cross within error bars at the critical temperature, giving a first indication that the 
percolation threshold coincides with the thermal critical point. The percolation probability 
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FIG. 3: Average Ising HT graph density as a function of the reduced inverse temperature (3/ (3 C 
for cubic lattices varying in linear size from L = 12 to L = 48. 

at the threshold we estimate to be Ps = 0.05(2). To obtain a more refined test, we apply 
finite-size scaling to this observable. When replotted function of {(3/(3 c - l)L lA \ with 
v the correlation length exponent of the 3D Ising model, the data is expected to collapse 
onto a universal curve. The inverse critical temperature and v have been determined to high 
precision in Refs. fl~2| IT3]. 

p c = 0.22165459(10), v = 0.63012(16), (20) 

respectively. Figure [5] shows that these values indeed produce a good collapse of the data 
in the entire temperature range. Since this is achieved without any adjustable parameter, 
we arrive at the important conclusion that the percolation threshold of the HT graphs 
coincides within error bars with the thermal critical point. Moreover, since the standard 
Ising correlation length exponent v has been used, it follows that the relevant diverging 
length scale for the HT graphs in the vicinity of the critical temperature is provided by the 
spin correlation length. 

We proceed to determine the fractal dimension D of the HT graphs. This can be done by 
measuring, for example, the percolation strength or the average graph size at the percolation 
threshold and applying finite-size scaling to the data obtained for different lattice sizes L. 
Since xg was found to show large corrections to scaling for small L, the observable P^, 
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FIG. 4: Percolation probability Ps of the Ising HT graphs as a function of the reduced inverse 
temperature /3/f3 c for cubic lattices varying in linear size from L = 12 to L = 48. 
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FIG. 5: Collapse of the data in Fig. [4] with the choices (20) for (5 C and v. 



showing only small corrections, is used to estimate D. Table [IT] summarizes the results 
of two-parameter fits using the nonlinear Marquardt-Levenberg algorithm for various fit 
intervals. On the basis of these fits, we estimate the critical exponent (3q/v to be 
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TABLE II: Percolation strength exponent 3gI v of the Ising HT graphs at the critical temperature 



(20) as obtained through two-parameter fits in the indicated intervals. 



L 


0c /u 


y 2 /DOF 

A / 


L 




Qclv 


y 2 /DOF 


8-64 


1 256(5) 


1.61 


10 


- 64 


1.265(7) 


1.35 


8 — 56 


1 955f'6) 


1 68 


1 


— 56 


1 964f7 N l 


1 43 

± .10 


8-48 


1.255(6) 


1.75 


10 


-48 


1.266(8) 


1.43 


8-40 


1.258(6) 


1.57 


10 


-40 


1.271(7) 


0.96 


8-32 


1.254(7) 


1.57 


10 


- 32 


1.267(8) 


1.05 


8-24 


1.253(7) 


1.50 


10 


- 24 


1 269(8) 


0.85 


8-20 


1.245(8) 


1.19 


10 


- 20 


1.260(11) 


0.83 


12 - 64 


1.256(5) 


1.61 


14 


-64 


1.270(10) 


1.41 


12 - 56 


1.270(9) 


1.41 


14 


-56 


1.269(12) 


1.53 


12-48 


1.273(9) 


1.37 


14 


-48 


1.273(13) 


1.51 


12-40 


1.281(7) 


0.68 


16 


-64 


1.261(13) 


1.37 


12 - 32 


1.279(9) 


0.79 


16 


-56 


1.258(15) 


1.47 


12 - 24 


1.285(7) 


0.34 


16 


-48 


1.261(17) 


1.49 



Pg/v = 1.2651(65), 



(21) 



corresponding to the largest possible fit interval (L = 10 — 64) that still gives a good fit 
quality with a x 2 P er degree of freedom (DOF), x 2 /DOF = 1.35. Figure [7] shows that, 
with this choice, the data collected in the vicinity of the critical temperature on lattices of 
different size L fall onto a universal curve when both axes are properly rescaled. 



The estimate (21) for the graph exponent leads to the estimate 

D = d-f3 G /u = 1.7349(65) 



(22) 



(d = 3) for the fractal dimension of the HT graphs at the critical point. This fractal 
dimension is closer to that of a self-avoiding walk in 3D, for which [14J D = \ jv = 1.7001(32), 
than to that of a Brownian random walk, for which D = 2. Given this estimate for the fractal 



dimension, the relation (17) leads to the estimate 

a = 0.9147(42) 
13 



(23) 




FIG. 6: Log-log plot of the percolation strength of the Ising HT graphs at the critical temper- 



ature (20 ) as a function of the lattice size L. The straight line is obtained through a two-parameter 



fit in the interval L = 10 — 64. 
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FIG. 7: Collapse of the Ising percolation strength data measured on lattices of different size L (left 
panel) when L^°' l 'P 00 is plotted as a function of (/3/f3 c - \)L X I V (right panel) on a semilogarithmic 
scale. 

for the graph distribution exponent a, which for self-avoiding random walks is unity. 
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IV. XY MODEL RESULTS 



We proceed by analyzing the HT graphs of the 3D XY model. The critical temperature 
of the truncated XY model, where links are allowed to be at most singly occupied, has to our 
knowledge not been determined before. To arrive at an accurate estimate of the percolation 
threshold, we consider the percolation strength data and search for the best data collapse 
given the value 

v = 0.6717(1) (24) 

for the XY correlation length exponent recently reported in Ref. [15J. This is done by 
rendering a motion picture out of about 300 single frames showing the data collapse for 
different values of K c . Successive frames correspond to slightly increased values of K c . A 
media player, such as MPlayer that can go forward and backward frame by frame, is used to 
play the motion picture, and to determine the value of K c with the best collapse. The quality 
of the collapse is established visually. Error estimates are based on the number of successive 
frames for which the quality of collapse remains roughly the same. We have checked this 
method by applying it to the Ising model, where the critical temperature is known to high 
precision, and obtained surprisingly good results. For the truncated XY model, we arrive in 
this way at the estimate 

K pcr = 0.22288(5), (25) 
which is to be compared to the value K c = tanh/3 c = 0.218095- •• of the Ising model. 



Figure |8| shows the collapse of the data achieved with the estimate ( 25 ) of the percolation 
threshold. 

As for the Ising model, we determine the fractal dimension of the HT graphs of the 



XY model by measuring the percolation strength P M at the percolation threshold (25) on 
lattices of different size. Table [TTT] summarizes the results of two-parameter fits using various 
fit intervals. On the basis of these fits, we estimate the critical exponent (3q/v to be 

P G /u = 1.2374(66), (26) 

corresponding to the largest possible fit interval (L = 10 — 64) that still gives a good fit 



quality (% /DOF = 1.21). Figure 10 shows that with this choice, the data collected in the 



vicinity of the percolation threshold on lattices of different size L fall onto a universal curve 



when both axes are properly rescaled. The result (26) leads to the estimate 
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FIG. 8: Collapse of the percolation probability Ps of XY HT graphs measured on lattices of 
different size L (left panel) when replotted as a function of (K/K pey - — l)L l / u (right panel). 



TABLE III: Percolation strength exponent (3q/v of the XY HT graphs at the percolation threshold 



(25) as obtained through two-parameter fits in the indicated intervals. 



L 




X 2 /DOF 


L 


Pg/v 


x 2 /dof 


6-64 


1.203(7) 


4.03 


8-64 


1.221(7) 


2.34 


6-56 


1.200(7) 


3.93 


8-56 


1.217(7) 


2.36 


6-48 


1.196(6) 


3.21 


8-48 


1.212(7) 


1.85 


6-40 


1.192(6) 


2.73 


8-40 


1.208(7) 


1.60 


6-32 


1.188(7) 


2.66 


8-32 


1.203(8) 


1.71 


6-24 


1.185(7) 


2.83 


8-24 


1.200(9) 


1.98 


6-20 


1.179(8) 


2.35 


8-20 


1.192(11) 


1.95 


L 


Pq/v 


x 2 /dof 


L 


Pg/v 


x 2 /dof 


10-64 


1.237(7) 


1.21 


12-64 


1.244(8) 


1.13 


10-56 


1.235(7) 


1.26 


12-56 


1.241(9) 


1.22 


10-48 


1.229(6) 


0.88 


12-48 


1.234(9) 


0.89 


10-40 


1.224(7) 


0.77 


14-64 


1.258(7) 


0.56 


10-32 


1.221(8) 


0.90 


14-56 


1.257(8) 


0.53 


10-24 


1.220(11) 


1.12 


16-64 


1.258(9) 


0.60 


10-20 


1.213(15) 


1.36 


16-56 


1.256(10) 


0.68 
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FIG. 9: Log-log plot of the percolation strength of the XY HT graphs at the percolation 



threshold (25) as a function of the lattice size L. The straight line is obtained through a two- 



parameter fit in the interval L = 10 — 64. 
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FIG. 10: Collapse of the XY percolation strength data measured on lattices of different size 
L (left panel) when L^^P^ is plotted as a function of (K/K per — l)L l / u (right panel) on a 
semilogarithmic scale. 



D = 1.7626(66) 



(27) 



for the fractal dimension of XY HT graphs at the percolation threshold. Comparison with 



the estimate (22 ) for the fractal dimension of the Ising HT graphs shows that the XY graphs 



are more crumpled although still much less so than a Brownian random walk. For the graph 
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distribution exponent a we arrive at the estimate 



a = 0.8446(45), (28) 



which is lower than the estimate (23) for the Ising model. Our estimate (27) is in good 
agreement with the value D = 1.7655(20) recently reported by Prokof'ev and Svistunov jl] 
for the |0| 4 theory which they obtained using their worm algorithm [6J. As stated in the 
Introduction, the HT graphs of that model as well as the worm update algorithm used to 
simulate them are completely different from the XY HT graphs and the plaquette update. 
Yet, despite these differences, the fractal dimensions of the two models, which share the 
same universality class, agree within error bars. We take this as a strong indication that, as 
expected, universality holds for the fractal structure of HT graphs. 



V. CONCLUSIONS 



In this paper, it was shown that the geometric Monte Carlo approach originally introduced 
in the context of 2D spin models [2], in which HT graphs are simulated directly and analyzed 
with the help of percolation observables, can also be applied to 3D spin models. The 3D 
Ising HT graphs were shown to percolate right at the critical temperature, which is known to 
high precision. The phase transition to the ordered, low-temperature state in this spin model 
was shown to be reflected by a proliferation of HT graphs. Also, through data collapse, it 
was shown that the diverging length scale relevant to the HT graphs in the vicinity of the 
percolation threshold is the spin correlation length. With the help of finite-size scaling, the 
fractal dimensions of the closed Ising and XY HT graphs were determined as in percolation 
theory. Both models are handled similarly, with the only difference that, in contrast to the 
Ising HT graphs, the XY HT graphs are oriented. Finally, it was shown that universality 
holds for the fractal structure of HT graphs. 
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